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ABSTRACT 

The distribution of orbital period ratios of adjacent planets in extra-solar planetary systems discov¬ 
ered by the Kepler space telescope exhibits a peak near ~ 1.5-2, a long tail of larger period ratios, 
and a steep drop-off in the number of systems with period ratios below ~ 1.5. We find from this 
data that the dimensionless orbital separations have an approximately log-normal distribution. Using 
Hill’s criterion for the dynamical stability of two planets, we find an upper bound on planet masses 
such that the most common planet mass does not exceed 10 _3 ' 2 m*, or about two-thirds Jupiter mass 
for solar mass stars. Assuming that the mass ratio and the dynamical separation (orbital spacings 
in units of mutual Hill radius) of adjacent planets are independent random variates, and adopting 
empirical distributions for these, we use Hill’s criterion in a statistical way to estimate the planet mass 
distribution function from the observed distribution of orbital separations. We find that the planet 
mass function is peaked in logarithm of mass, with a peak value and standard deviation of log m/M® 
of ~ (0.6 — 1.0) and ~ (1.1 — 1.2), respectively. 


1. INTRODUCTION 

What is the mass distribution function of plan¬ 
ets in the Universe? In the past two decades, 
more than 1700 exoplanets have been discovered; of 
these, 502 have measured masses (Exoplanet Archive, 
http://exoplanetarchive.ipac.caltech.edu/, as of Septem¬ 
ber 16, 2014). The distribution of these planet masses 
is shown in Fig. [T] We observe two local peaks in this 
apparent distribution, near ~ 1 Mj and near ~ 10 M®: 
Jupiter-mass and “super” Earth-mass planets are com¬ 
mon in the discovered population of planets. How¬ 
ever, this apparent mass distribution suffers from many 
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Fig. 1.— The distribution of log-mass of con¬ 
firmed exoplanets with measured masses (data from 
(http://exoplanetarchive.ipac.caltech.edu/, retrieved on Septem¬ 
ber lOj, 2U14J. t he black points indicate the masses of the solar 
system planets. Note that this is a semi-log plot. 
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difficult-to-quantify selection biases, so we must exercise 
great caution in interpreting these features. At the very 
least, they are of unknown significance for the intrinsic 
distribution of planet masses. 

In some contrast with the less-than-30% fraction 
of all exoplanets for which we have measured masses, 
the orbital periods of nearly 100% of exoplanets are 
quite well determined - indeed, the periodicity of the 
orbital motion of planets is predominantly how they are 
discovered. The Kepler miss ion, current ly the la rgest 
systematic exoplanet survey dBorucki et al.1 IJOTTI) . has 
provided a wealth of data on planets and planetary 
systems in the Galaxy. A large subset, about 65%, of all 
confirmed exoplanets are found in planetary systems har¬ 
boring two or more planets. Significantly, several studies 
of the Kepler data on multiple-planet systems have 
concluded that planetary systems are coplanar to within 
a few degrees (iLissauer et all 1 1201 It iTremainc fc Pond 


2012t iFang fc Margod 20121: Johansen et all l2012t 
Figueira et all 120121 iFabrvckv et al.1 2014), and that 


they are likely closely-packed dFang fc Maraotll2013l) . 

Here we leverage the Kepler data on orbital periods in 
multiple planet systems, together with theoretical under¬ 
standing of the long term dynamical stability of coplanar 
planetary systems, to estimate the planet mass distri¬ 
bution function. First, we use the observational data 
on orbital periods to compute the distribution of orbital 
spacings in systems of N > 2 planets. We then reason 
that it must be possible to deduce planetary masses from 
the observational data of orbital period ratios if orbital 
spacings are determined by long term dynamical stabil¬ 
ity. Numerical studies have shown that the relationship 
between orbital spacing and planet masses is necessar¬ 
ily statistical in nature because multiple planet systems 
exhibit chaotic dynamics. We adopt the ansatz that the 
orbital spacing measured in units of the mutual Hill ra¬ 
dius —the so-called “dynamical separation”— is a ran¬ 
dom variate, and we adopt an empirical distribution for 
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this parameter. This leads us to estimates of the dis¬ 
tribution of the total mass of adjacent planets relative 
to the stellar host mass. We then consider two limiting 
cases, that adjacent planets have a random mass ratio 
or that they tend to be similar to each other in mass. 
Finally, we convolve with the observed distribution of 
stellar host masses to convert from planet-to-star mass 
ratios to planet masses, to calculate the planet mass dis¬ 
tribution function. 

Our approach assumes that dynamical separations, to¬ 
tal mass and mass ratios of adjacent planets are inde¬ 
pendent random variates, and we neglect any corelations 
with stellar host mass and the age of the system. With 
these simplifications, we arrive at a theoretical estimate 
of the planet mass function based only on the observa¬ 
tional data of the orbital periods of exoplanets and the 
masses of their stellar hosts. We also make indepen¬ 
dent estimates of the planet mass function based on em¬ 
pirical mass-radius relations and observational data of 
planet radii, and we compare these with the dynamical 
stability-based estimate. 

The true mass distribution function will of course be in¬ 
creasingly better determined with ongoing observational 
efforts to measure the masses of a large population of 
extra-solar planets by means of complementary tech¬ 
niques (transits, radial velocities, astrometry, etc.). Our 
simple theoretical prediction for the planet mass distri¬ 
bution function may be useful for the interpretation of 
such forthcoming observations. Our work may also be 
useful for the planning and interpretation of numerical 
studies of the dynamical stability of planetary systems. 
Knowledge of the planet mass function is important for 
understanding the physics of planet formation in differ¬ 
ent mass regimes, as well as for assessing the abundance 
of planets like our own home planet. 


2. ORBITAL SPACINGS 


We will use the following notation: m* is the stellar mass; 
the planets’ masses, orbital semimajor axes and orbital 
periods are rrii,ai and Ti, with i = 1 , 2 ,... in order of 
increasing distance from the host star. We define the 
mass ratios, 



A H = 


rrij + m i+ i 
m* 


min{mj,m i+ i} 
max{m,, mq i} 


(1) 

( 2 ) 

( 3 ) 


Note that 0 < ^ < /ij C 1 and 0 < 7 * < 1. The 
masses of the two nearest-neighbor planets are given by 

7i(l + 7i) -1 Ai TO * an d (1+7 

We also define the period ratio for nearest-neighbor 
planets, Vi = T i+ i/Ti. In Figure [2j we plot the dis¬ 
tribution of Vi for the ensemble of 373 multiple planet 
systems discovered by Kepler, there are 566 period ra¬ 
tios in this data. It is a broad distribution, with paucity 
of period ratios close to 1 , a peak near 1 . 6 , and a long 
tail of large values. There is also interesting fine struc¬ 
ture within this broad distribution, particularly a trough- 
peak feature near low order resonant period ratios, such 
as 3/2 and 2/1, that has been d i scussed in several re- 
cent papers dLithwick fc Wull2012t iBatvein fe Morbidellil 



Fig. 2.— The period ratio distribu tion in multiple planet systems 
discovered by Kepler. (Data from [Fabrvckv et ali (120141) .) The 
dot-dashed curve is a smoothed version of the histogram (smoothed 
with a Gaussian kernel). The black points indicate solar system 
values. The vertical dotted lines indicate locations of low order 
resonant values (3/2, 5/3, 2/1,7/3, 8/3). 


120131 iPctrovich ct al.ll2013l iFabrvckv et al.ll2014lf . In the 
present work, we attempt to understand the overall dis¬ 
tribution of V . 

Examining Fig. [2] it is not difficult to be persuaded 
that the steep drop in the V distribution from near 
P « 1.5 to P « 1.3 and the paucity of systems with 
V close to unity is likely owed to the instability of very 
closely spaced orbits, due to mutual planetary perturba¬ 
tions: larger planet masses require larger orbital spacings 
for long term dynamical stability; higher planet multi¬ 
plicity and higher orbital eccentricities also would tend 
to require larger orbital spacing for stability (with the 
exception of librating resonant orbits). Therefore, if the 
orbital spacings in multi-planet systems are related to 
their long term dynamical stability, the planet masses 
must be related to the orbital spacings. By use of Ke¬ 
pler’s third law, the dimensionless orbital spacing, T>i, is 
related to the period ratio of adjacent planets, 


Vi = 2 


Qi+l — 
Oi+1 + di 


V? - 1 
V? + 1 


( 4 ) 


Note that 0 < Vi < 2. 

The observed distribution of log V for Kepler planets 
is shown in Figure [3] We find that a Gaussian function, 
with mean xt> = (log V) = —0.318 and standard devia¬ 
tion a-D = 0.231, fits the data fairly well; a y 2 test gives 
p-values > 0.1. Of course, the Gaussian must be trun¬ 
cated at a maximum value, log V = log 2. Therefore, 
formally, the best-fit probability density function (PDF) 
for X = log V is expressed as 


F-d(x) 


_i_cxnf 

P[ 2<r£ J 

0 


if x < log 2 
otherwise, 


( 5 ) 


where 


$23 


j. ( log2 -^), 

(723 


( 6 ) 


and $(•) is the cumulative distribution function of the 
standard normal distribution. For the computed values 
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Fig. 3.— The distribution of the orbital spacing, X> (Eq. QJ, of 
adjacent planets in multiple planet systems discovered by Kepler. 
The dot-dashed curve is the best-fit Gaussian function. The black 
points indicate solar system values. 

of the mean and standard deviation, we find $£> ~ 0.996. 
Therefore we incur only a very small error in approximat¬ 
ing Fx>(x) as an untruncated Gaussian function. 

In other words, the PDF of the dimensionless orbital 
spacing, T>. is a nearly log-normal distribution. A log¬ 
normal is a skewed distribution; the mean, median and 
mode of V are 0.554, 0.481 and 0.362, respectively. Note 
that a log-normal distribution can resemble a power-law 
distribution over a fairly wide range of the parameter 
away from the peak. 

We remark that the dimensionless orbital separations 
of solar system planets (indicated by the black dots in 
Fig. [3j) are not dissimilar to those of the Kepler planets. 

3. DYNAMICAL STABILITY 

Let’s consider the dynamical stability of a pair of ad¬ 
jacent planets in the simplest case, that of a two-planet 
system. For nearly co-planar and initially circular orbits, 
the smallest orbital spacing that is dynamically stable is 
given by Hill’s criterion: 

V = K^y, with K = 2>/3. (7) 

Note that this stability criterion is insensitive to how the 
total planet mass is partitioned between the two planets. 
It is independent of the actual distance of the planets 
from the host star. 

For systems with more than two planets, there is no 
known analytic criterion for dynamical stability, but we 


TABLE 1 

Solar system planets 


Planet pair 

A 

V 

7 

K 

Venus- Mercury 

2.61E-06 

2.55 

0.678E-01 

63.4 

Earth-Venus 

5.45E-06 

1.63 

0.815 

26.3 

Mars-Earth 

3.33E-06 

1.88 

0.107 

40.1 

Jupiter-Mars 

9.55E-04 

6.31 

0.338E-03 

16.0 

Saturn-Jupiter 

1.24E-03 

2.48 

0.299 

7.90 

Uranus-Saturn 

3.29E-04 

2.85 

0.153 

14.0 

N ept une- U ranus 

9.52E-05 

1.96 

0.848 

14.0 


expect that in close-packed systems of N > 3 planets, 
orbital spacings must exceed those required by Hill’s cri¬ 
terion, i.e., K must exceed 2-\/3. We can look to nu¬ 
merical studies of dynamical stability of planetary sys¬ 
tems for insights. Several such investigations have been 
published, many focussed on particular systems, but a 
few on the broader theoretical question of the dynami¬ 
cal lifetimes of multiple planet systems as a function of 
planet masses and orbital spacings. Numerical results for 
nearly coplanar, low eccentricity multi-planet systems (of 
3 < IV < 20 equal-mass planets with 10 -9 < /.q < 10 -5 ) 
show that dynamical stability times in excess of ~ 10 s Ti 
require that adjacent planet pairs must have K > 8 
and that this lower limit on K is only weakly dependent 
on planet mass and planet m ultiplicity dChambers et al.l 
119961 ISmith fe Lissauerll2(j09l l. For higher planet masses 
(10 -3 ' 4 < pi < 10 ~ 2 ' 4 ), the critical K is some¬ 
what smaller, K > 5 dMarzari fe Weidenschillind 120021 
IChatteriee et all[200811 . The notation “A”, “/?” and K 
has been used by various authors for the same parame¬ 
ter; here we have adopted K. This parameter, which is 
the orbital separation in units of the mutual Hill radius, 
is often called “dynamical separation”. 

An important insight from these studies is that the 
chaotic nature of multiple planet systems necessitates a 
statistical description of the relationship between planet 
masses, orbital separations and the system’s dynami¬ 
cal stability time, i.e., that K depends sensitively on 
initial conditions and is better described as a random 
variate. The distribution of K depends upon the plan¬ 
etary architecture. For the singular but well-studied 
example of the solar system, investigations of its long 
term dynamics have concluded that it is marginally sta¬ 
ble on timescales comparable t o its age (iLaskad 1996; 
Haves! 120081 iHaves et al.1 120101 : iLithwick fc Wul 2014 : 
Batygin et al.ll2015h . Taking adjacent pairs of solar sys¬ 
tem planets, we find that K ranges from ~ 8 (for Jupiter- 
Saturn) to ~ 63 (for Mercury-Venus), with a mean value 
of 26 (see Table 1). These values are larger than required 
by the two-planet Hill’s stability criterion, and reflect the 
effects of non-trivial orbital eccentricities, mutual incli¬ 
nations, unequal planet masses and planet multiplicity. 
With this as motivation, we adopt a heuristic criterion 
for the long term stability of systems with more than 
two planets: a straightforward generalisation of Eq. [7] in 
which we treat K as an independent random variate. We 
will denote by Pk {-) the PDF of AT, and denote by Fk(-) 
the PDF of log K. 

4. ANALYTICAL ESTIMATES OF THE DISTRIBUTION OF 
TOTAL MASS OF ADJACENT PLANET PAIRS 

The following rearrangement of Eq.[7]will be useful for 
our analysis, 

log fi = 3(log V — log K) + log 3. (8) 

We can then estimate the probability density function 
Fp(-) of log /i, as a convolution of the PDFs of log V and 
of log AT, 

F -( x )= f-oo dx vJZo dx KF K (x K )F D (xv) (Q) 
M xS(x — 3x-d + 3xk — log 3), 

where S(x) is the Dirac delta function. One challenge 
we face is that the published numerical studies have 
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not reported the PDF of AT, nor have they explored 
the large parameter space of systems containing unequal 
mass planets. This requires systematic numerical inves¬ 
tigation, which is feasible with modern computers, but 
has yet to be undertaken. In the absence of such knowl¬ 
edge, we will adopt some plausible ansatzs for the PDF 
of AT. 

As a simple illustration, let’s first consider Fk(xk) = 
8{xk — log AT*). Then, with the best-fit Gaussian func¬ 
tion for the PDF of log V (Eq.0, it is straightforward to 
determine that the distribution of log jl is also Gaussian, 
with mean (log jl) = —0.48 — 3 log AT*, and standard de¬ 
viation ap = 0.693. If we set AT* = 2-\/3, the minimum 
value needed for dynamical stability of two planets, then 
log jl has a Gaussian distribution of mean (log jl) = —2.09 
and standard deviation 0.693. The corresponding dis¬ 
tribution of jl is log-normal, with median value 10 -2 09 
and mode 10 -3 ' 20 . Because this describes the maximal 
jl that is dynamically stable, we can conclude that the 
most common planet mass does not exceed 10 _3 20 to», 
or about two-thirds Jupiter mass for solar mass stars. 

For a more realistic estimate, let’s consider a broad 
distribution of K values. We are inspired to consider 
a Gaussian distribution for its simplicity, because it fol¬ 
lows from Eq. [8] that log jl has a Gaussian distribution if 
log V and \ogK are both Gaussian variates. We are also 
motivated by the AT values in the solar system (Table 1): 
log AT is in the range 0.9 to 1.8, with mean (log AT) = 1.32 
and standard deviation ctk =0.31. Although the sample 
is small, its log AT distribution is not inconsistent with a 
normal distribution. For Fk(xk ), we therefore adopt a 
Gaussian functiorQ with the mean and standard devia¬ 
tion chosen to match the values found for the solar sys¬ 
tem. Then, it is straightforward to compute that the re¬ 
sulting Gaussian PDF for log jl has mean (log jl) = —4.44 
and standard deviation ap = 1.16. The median value of 
jl in this case is 10“ 4 ' 44 ; the mode is ~ 10~ 7 ' 5 , about 
three orders of magnitude smaller than the upper bound 
derived above with AT* = 2y/Z. 

It is evident that dynamical stability implies that the 
PDF of log jl has small probability density at both small 
and large values and a peak at an intermediate value. 
This shape is inherited from the nearly Gaussian distri¬ 
bution of log 2?, which in turn is derived from the ob¬ 
served distribution of period ratios. In particular, the 
drop-off at small masses is inherited from the steep drop¬ 
off in the number of observed systems with period ratios 
smaller than ~ 1.5. 

We remark that the mean value of log jl decreases with 
increasing mean K. However, the peaked shape of the 
PDF of log jl is not severely dependent on the particular 
functional form of the K distribution that we adopted. 
For example, a flat, uniform random distribution of AT in 
a range K min > 2y/3 to K max = 95 (slightly wider than 
the range found in the solar system) also yields a peaked 
distribution of log jl. The skewness of the jl distribution 
does depend on the dispersion of the K distribution, and, 

1 Formally, the PDF of log A' must vanish for K < 2\/3, to 
satisfy the Hill criterion (Eq. [FJi for the limiting case of two planet 
systems. This means that we should adopt a truncated Gaussian 
PDF, analgous to the case of the PDF of logD (Eq(SJ. For the 
parameters of interest here, the normalization factor is nearly unity, 
so we neglect the truncation. 


consequently, the most probable value of jl depends upon 
this dispersion as well. Our choice of the PDF of K is of 
course motivated by the solar system and could be con¬ 
sidered biased; we discuss this point further in Section 6 . 

5. NUMERICAL ESTIMATES OF THE PLANET MASS 
DISTRIBUTION FUNCTION 

To determine the distribution of individual planet 
masses requires additional considerations: are the masses 
of adjacent planet pairs correlated or are they indepen¬ 
dent? 

Let P^,{x) be the probability density function of /x; 
(planet mass as fraction of stellar mass). Let’s consider 
the limiting case in which adjacent pairs of planets have 
mass ratio which is a fixed constant, 7 j = 7 *. Then the 
PDF of /i,; is straightforwardly derived from that of jl, 

p n(x) = r[( 1 +7*) p a(( 1 +7*)2;)+( 1 +7* _1 )^/i(( 1 +7*“ 1 )^)]- 

( 10 ) 

On the other extreme, if the masses of adjacent planet 
pairs are considered independent, then we have the fol¬ 
lowing relationship between P ^ and Pp, 

p p( x )=[ dy P»(y)Pn(x-y). (11) 

J 0 

We will not attempt to solve this implicit equation for P ^, 
but we can note that, depending on its functional form, 
the mass ratio, 7 , of planet pairs may or may not be 
independent of jl, even if the masses of adjacent planets 
were independent. 

In reality, adjacent planet masses are likely to be nei¬ 
ther perfectly correlated nor perfectly independent. We 
therefore consider 7 to be a random variate, and for sim¬ 
plicity, we assume that it is independent of the total mass 
of the planet pairs. We consider two cases: (i) a uniform 
PDF for 7 , P-y(x) = 1 or (ii) a PDF with peak at 1. For 
the latter, we chose P 7 to be a Gaussian with mean 1, 
standard deviation 0.3, truncated at 0 on the left and 1 
on the right; in this distribution, approximately half of 
all planet pairs have 0.8 < 7 < 1). This distribution is 
similar to the assumption that neighboring planets tend 
to have similar masses, whereas the uniform distribution 
allows any value of their mass ratio with equal probabil- 
ity. 

We carry out the numerical calculation of the planet 
mass distribution function as follows. Starting with the 
observational data of the period ratios, Vi, of the Kepler 
planets, we first calculate Vi. Then we calculate jii with 
the help of Eq. [ 8 ] and a random value of K from its pre¬ 
scribed PDF; we adopt a Gaussian PDF of log A' with 
mean 1.32 and standard deviation 0.31, as discussed in 
the previous section. Next we compute the individual 
Pi s with the help of Eq. [2] [3] and a random value of 7 
from its prescribed PDF. We average over 1000 realiza¬ 
tions of the random choices of log AT and 7 . We take 
one additional step: to compute the individual planet 
masses, rrii = ptm*, we adopt the stellar host masses 
of the Kepler multiple-planet systems (obtained from es- 
ti mates of the st e llar su rface gravity an d stellar radius 
b y iBatalha et al.l (1201311 . as reported in iFabrvckv et al.l 

{MIT 

The results are shown in Fig. [4], where we plot in con¬ 
tinuous line the case of 7 uniform random distribution 
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Fig. 4. — The distribution of log-mass of Kepler planets, derived 
from their period ratios and a heuristic criterion for dynamical sta¬ 
bility. The continuous line curve and dot-dashed curve are the re¬ 
sults obtained by assuming a uniform random and a half-Gaussian 
distribution, respectively, of the ratio of adjacent planet masses, 7 
(Eq.[3)). The black curves are based on the observed period ratios, 
while the blue curves are based on the de-biased distribution of 
period ratios. The black points indicate the masses of solar system 
planets. 

on ( 0 , 1 ), and in dot-dashed line the case of 7 having a 
half-Gaussian PDF peaked at 1. We observe that (a) 
the estimated planet mass function is not very sensitive 
to the choice of P 7 , and (b) the PDF of the logarithm 
of planet mass does not increase monotonically as the 
mass decreases. The distribution of logm/M® is found 
to be peaked, with mean 0.64(0.72) and standard de¬ 
viation 1.21(1.17) for 7 uniform random (half-Gaussian 
peaked at 1 ). 

The above analysis is based on the observed orbital 
period ratios of Kepler planets. We can ask about the 
selection biases and incompleteness of the observed dis¬ 
tribution of V. If a higher sensitivity survey, analogous 
to the Kepler survey, were to be carried out, we might ex¬ 
pect that smaller planets would be discovered in greater 
numbers, over the same range of orbital periods (limited 
by the length and cadence of the survey). It is not im- 
medi ately obvious h ow this would affect the distribution 
of V. iSteffenl (120131) has compared the distribution of V 
of the high multiplicity and the low multiplicity systems 
within the Kepler sample and concluded that the period 
ratio distribution of Kepler’s harvest of multi-planet sys¬ 
tems is a fair sa mple of the intrinsi c distr ibution, at least 
for V < 5 or 6 . iSteffen fe Hwand (120151) have analyzed 
the incompleteness of the observed V distribution due to 
planets missed by the Kepler data reduction pipeline as 
well as due to geometric bias against detection of non- 
coplanar planets. The authors assumed that the mutual 
inclinations of planetary orbits have a Rayleigh distri¬ 
bution with width parameter a = 1.5 degrees. Their 
debiased V distribution is broadly similar to the ob¬ 
served distribution, but has higher probability densities 
for larger period ratios (see their Figure 4). Using this 
debiased V distribution, we repeated our calculation of 
the planet mass distribution; the results are plotted in 
blue in Figure |U We find that the resulting distribu¬ 
tion of logm/M® has mean 0.91(0.98) and standard de¬ 
viation 1.18(1.14) for 7 uniform random (half-Gaussian 
peaked at l). Steffen (2015, personal communication) 


also provided us with a debiased V distribution based 
on a model Rayleigh distribution of the mutual inclina¬ 
tions with width parameter a = 3°; the resulting planet 
mass distribution is insignificantly different than that for 
a = 1.5 degrees. 

In summary, we find that the planet mass distribu¬ 
tion derived from use of the heuristic dynamical stability 
criterion and the observed and debiased distribution of 
period ratios of Kepler planets is peaked in logarithm of 
planet mass, with a peak value of log m/M® of 0.6 — 1.0, 
and standard deviation 1.1 — 1.2. 

6 . COMPARISON WITH OTHER ESTIMATES 

All the planets discovered by Kepler have fairly well de¬ 
termined values of their planetary rad ii, with uncertain¬ 
ties of about 30% dSilburt et all 120151 ). Several studies 
have proposed em piric al relations betwee n planet masses 
and planet radii. iLissauer et ah : (|2011l) have proposed 
the following, based on the well-known properties of the 
solar system planets: 

m/M® = (R/R®)“, (12) 

where M® and R® are the mass and radi us of Earth, a = 
2.06 i f R > R® and a = 3 if R < R®. IWu & Lithwickl 
( 2013 ) have proposed a slightly different relation, based 
on observational data of a subset of Kepler planets whose 
masses have also been determined observationally: 

m/M® = 3(R/R®). (13) 

IWeiss fc Marcvl (120141 ) report a slightly different empir¬ 
ical mass-radius-density relation, based on an updated 
list of Kepler planets of radius R < 4R® whose masses 
have been measured; this can be expressed as follows: 

(0.441 + 0.615R/R®)(R/R®) 3 if R < 1.5R®, 
2.69(R/R®) 0 ' 93 if 1.5R® < R < 4R®. 

(14) 

Using each of these mass-radi us relations and t he da ta 
for planetary radii (reported in iFabrvckv et ah! (j2014f )). 
we computed the masses of the 939 Kepler planets in 
multiple-planet systems, and then used Gaussian kernel 
density estimation to compute the PDF of the logarithm 
of masses. For R > 4R®, we supplemented Eg. ITT! with 
Eq.HH 

As an aside, we note that, having computed the planet 
masses by using mass-radius relations, we then computed 
K values for the adjacent planet pairs. From these val¬ 
ues, we estimated the PDFs of log AT using Gaussian ker¬ 
nel estimation; these are shown in Fig. [5] We find that 
the three empirical mass-radius relations yield similarly 
peaked distributions of log AT, with mean values 1.32, 
1.29 and 1.31, and standard deviations 0.24, 0.23 and 
0.23, respectively. The mean values are similar to that 
of the solar system planets; the standard deviations are 
somewhat smaller. This provides support for the PDF 
of log K that we adopted by ansatz. 

The mass distributions obtained by use of the mass- 
radius relations (Eq. [T2l Eq. [131 Eq. [14]) are shown in 
Fig. [6] We see that these yield strongly peaked PDFs of 
logm/M®, with mean values 0.66, 0.80 and 0.68, respec¬ 
tively. For comparison, we also plot our theoretical esti¬ 
mates based on dynamical stability. The PDFs derived 
independently from the mass-radius relations peak near 
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log K 


Fig. 5. — Comparisons of the distribution of log K , where K is the 
orbital spacing in units of the mutual Hill radius. The dark grey, 
light grey and dot-dashed cu rves are the results o btai ned from the 
mass-radius relations of Eq. 1121 Eq. Il3l and Eq. 1141 respectively. 
The black continuous line curve is the Gaussian distribution that 
we adopted, with mean and standard deviation matching that of 
the solar system planets. The black points indicate solar system 
values. 



log m/M e 

Fig. 6 .— Comparisons of PDFs of the logarithm of planet mass: 
The dark grey, light grey and dot-da shed cu rves are the res ults 
from the mass-radius relations of Eq. Il2l Eq. Il3l and Eq. I 111 re¬ 
spectively. The black and blue continuous and dotted curves are 
our theoretical estimates (as in Fig. 0 . 

similar values to our dynamical stability-based PDFs. 
The latter have significantly larger dispersion, however. 
The smaller dispersions of the former are not entirely 
surprising, as the mass-radius relations describe empiri¬ 
cal best-fits and do not reflect th e uncertainties and dis - 
persion in the observational data. IWeiss fc Marcvl (I2014T ) 
note that the observational data have significant scatter 
about their empirical mass-radius relation, and that the 
scatter is not merely due to observational errors but may 
reflect intrinsic compositional diversity of planets. 

Overall, in comparison with the estimates based on 
mass-radius relations, the dynamical stability-based es¬ 
timate has a shallower power law slope of the planet mass 


function at planet masses larger than about ten Earth- 
mass and a steeper slope for masses below about one 
Earth-mass. These discrepancies may reflect an over¬ 
estimate of the width of the empirical distribution of 
log K adopted in our dynamical stability-based estimate, 
or they may reflect observational incompleteness (due to 
undetected small and/or long period planets) and the 
intrinsic scatter about the best-fit mass-radius relations. 
Future studies can test these alternate hypotheses. 


7. DISCUSSION AND SUMMARY 

A quantitative description of the architectures of plan¬ 
etary systems requires at least the following: the de¬ 
gree of orbital co-planarity, the distribution of their or¬ 
bital periods and spacings, and the distribution of plan¬ 
etary masses. The large number of planetary systems 
discovered by the Kepler mission allows a statistical 
assessment of these properties that have long eluded 
the theory of the formation and evolution of planetary 
systems. Several studies of the Kepler data suggest 
that planetary systems are flat, to within a few de¬ 
grees , similar to our ow n solar system dLissauer et all 
20T1 iTremaine fc Don til 120121: iFang fc Margotl 120121: 
Fabrvckv et al.T2014f l. Orbital sp acings have als o been 


the subject o f several studi e s dLissauer et al.l 120111: 
Steffenl (20131 : I Fabrvckv et al.1 12014 : Steffen fc Hwand 


20151/ as hav e planet masses dHoward et al.l 120101 : 
Wu fc Lithwickl 1 201. 'll IWeiss fc Marcvl I2014H . In the 


present work, we studied the distribution of orbital spac¬ 
ings and used dynamical stability to estimate planet 
masses in a statistical way. 

Regarding their orbital spacings, we found that the Ke¬ 
pler planets’ dimensionless orbital spacings, T> (Eq. [4]), 
have a nearly log-normal distribution (Fig. [3|. How 
can we understand this distribution of T>! In qualita¬ 
tive terms, a log-normal distribution is generated in a 
multiplicative random process, that is, when a positive 
definite variable v suffers random increments in propor¬ 
tion to its value, Sv oc v. If the successive increments are 
independent and large in number, then by the central 
limit theorem log v will be approximately normally dis¬ 
tributed. We conjecture that the log-normal distribution 
of V arises in the late stages of the dynamical evolution 
of planetary systems when the planetar y architecture is 
shape d (or re-shaped) by secular chaos dLithwick fc Wul 
1201411 . or possibly it arises in an earlier stage when a small 
number of planets emerge from their natal protoplane¬ 
tary disk but still embedded in a leftover disk of a large 
number of planetesimals. As a consequence of mergers or 
ejections of planetesimals and/or planets, the surviving 
planets undergo a random walk of their orbits; unsta¬ 
ble configurations are steadily winnowed. Our u nder¬ 
standing of this evolution is st i ll at an early stage |Rein| 
20121: lHansen fc Murravl 1 201. 'll iMinton fc Levisonl 12014 


Hands et al.l 12014 iChatteriee fe Fordl 1201511 . In future 
studies, it would be very useful to examine quantitatively 
how (or if) secular chaos and/or planetesimal-aided or¬ 
bital migration lead to a log-normal distribution of or¬ 
bital spacings. 

Regarding the masses of planets, we reason that it must 
be possible to deduce planetary masses from the obser¬ 
vational data of orbital period ratios if orbital spacings 
are determined by dynamical stability. However, there is 
no direct way to do so, therefore we made several simpli- 
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fying assumptions and ansatzs. First, we used the two- 
planet Hill’s stability criterion to derive an upper bound 
for the most common planet mass. Then, we general¬ 
ized Hill’s criterion in a statistical way for the stability 
of multi-planet systems (Eq. [ 8 } to compute the planet 
mass function. We assumed that the dynamical separa¬ 
tion (i.e., the orbital separation in units of the mutual 
Hill radius) and the mass ratio of adjacent planets are 
both independent random variates. We adopted plausi¬ 
ble distribution functions for these two parameters, based 
on our understanding of solar system dynamics. These 
empirical distributions and the assumed independence of 
the variables are admittedly major simplifications. These 
simplifications can be relaxed in a future study by de¬ 
termining the joint probability density function of K , ft 
and 7 by means of large scale numerical simulations of 
the dynamical stability of multiple planet systems. 

The multi-planet systems discovered by Kepler have 
been described as being rather unlike the Solar system, 
because their orbital periods are ~ 10 days and their 
masses are on the order of a few Earth masses; such 
planets are absent in the solar system. However, the 
dimensionless orbital spacings and the dynamical sepa¬ 
rations of solar system planets are not dissimilar to those 
of the Kepler systems (Figs. [3] and [5]), when we compute 
the latter independently from the observational data for 
planetary radii by using mass-radius relations. By these 
scaled measures of “planetary system architecture”, the 
solar system does not appear to be an outlier. 

A deeper look at the dynamical separations in the so¬ 
lar system shows that the terrestrial planets, Mercury- 
Mars, have K values significantly larger than those of 
the outer planets, J upiter-Neptune (s ee Ta ble 1). This 
was pointed out in llto fc Tanikawal (|2002h . It is also 
notable that the K values for the giant planets (Jupiter- 
Neptune) are closer to the numerically determined min¬ 
imum value, AT* ss 5 — 8 , that is necessary for the dy¬ 
namical stability of TV > 3 equal mass planetary systems 
for timespans of the age of the solar system, whereas 
the K values for the terrestrial planets are significantl y 
higher. We conjecture, following llto fc Tanikawal (12002 ). 
that this dichotomy is owed to the property that the 
subsystem of the four terrestrial planets is subject to 
long term dynamical excitation by the giant planets. 
T his is support ed by the numerical experiments reported 
in lHaves et al. j (2010) who found that the terrestrial plan¬ 
ets exhibit a much lower level of long term chaos if the 
gravitational perturbations of the giant planets were ab¬ 
sent. It is also possible that K is related to “dynamical 
age”, i.e., the age of the system measured in units of the 
orbital period of the innermost planet. We observe that 
in the solar system, the terrestrial planets’ dynamical 
age is 1 to 2 orders of magnitude larger than that of the 
outer planets (since the orbital period of Jupiter is about 
50 times that of Mercury). If so, then solar system-like 
planetary architectures may be better modelled with a 
bimodal PDF of K. For example, we could consider a 
PDF of log K consisting of an equally-weighted sum of 
two Gaussian functions having mean (log K) of 1.06 and 
1.61, respectively, and each having standard deviation 
0.2; solar system K values of the giant planets and ter¬ 
restrial planets, respectively, are roughly consistent with 
these parameters. With this choice of bimodal PDF of 
log AT, we repeated the numerical calculations described 


in Section 5 to compute the associated planet mass func¬ 
tion. The resulting distribution of log m is bimodal. The 
two peaks are near log m/M® values of —0.2 and 1.4; 
the overall mean value of log m/M® is 0.60(0.67) and 
the standard deviation is 1.41(1.37) for 7 uniform ran¬ 
dom (half-gaussian). The two peaks are of only slightly 
differing heights, and this PDF can also be described as 
approximately a plateau in the range —0.2 to 1.4. We 
do not belabor the results of this numerical experiment, 
however, because at this level of detail, we must also 
pay attention to possible correlations amongst the pa¬ 
rameters, K , ft, 7 and orbital periods. A comprehensive 
numerical study of the dynamical stability of multiple- 
planet systems can provide an improved estimate of the 
joint probability distribution of these parameters. 

In any study of the mass distribution of planets, the 
question arises of the definition of “planet”, a question 
that has been debated in both public and scientific fo¬ 
rums. At the upper end of the planet mass range, the lit¬ 
erature makes a distinction between “giant” planets and 
“brown dwarfs” near a mass of about 13Mj(« 4134M®). 
In the solar system, “planet” masses range from 0.055M® 
to 318M®, and, at the lower end of the mass range, the 
literature also recognizes “dwarf planets”, the most mas¬ 
sive being ~ 0.002M®. The results of our dynamical 
stability-based estimate of the mass function are most 
pertinent for the mass range spanning a few percent of 
an Earth-mass to a few hundred Earth-masses. Although 
these results can be smoothly extended beyond these lim¬ 
its, it is likely that different physical processes shape the 
mass function near these upper and lower bounds. 

Our results and conclusions are summarized as follows. 

1. The observed period ratios of adjacent planet pairs 
in multiple-planet systems discovered by Kepler in¬ 
dicate that the distribution of the dimensionless 
orbital separation (Eq. [4|) is approximately a log¬ 
normal function. 

2. The minimum dynamical separation (the orbital 
separation of adjacent planets in units of their mu¬ 
tual Hill radius, Eq. 0 . K = 2^, necessary for 
long term dynamical stability of two-planet sys¬ 
tems implies that the most common planet mass 
does not exceed 10 _ 32 m*. For solar mass stars, 
this is about two-thirds the mass of Jupiter. 

3. For plausible distributions of K and of the mass 
ratio, 7 (Eq. [3j, of adjacent planet pairs, our the¬ 
oretical estimate of the planet mass distribution 
function is peaked in log to. It is only weakly sen¬ 
sitive to the distribution of 7 . We estimate that the 
most probable value of logm/M® is ~ (0.6 — 1.0), 
and the standard deviation of the distribution of 
log m/M® is about 1.2. 

4. The planet mass distribution computed indepen¬ 
dently from the observational data on planetary 
radii (by use of empirical mass-radius relations) is 
also peaked in logm at similar peak values, but 
has smaller dispersions. These discrepancies may 
reflect an overestimate of the width of the distribu¬ 
tion of log K adopted in our theoretical estimate, 
or they may reflect observational incompleteness of 
the measured masses at low and high planet masses 














and the intrinsic scatter of the masses about the 
best-fit mass-radius relations. Future studies can 
test these alternate hypotheses. 

5. In deriving the dynamical stability-based esti¬ 
mates, we assumed that K and 7 are independent 
random variates. We adopted PDFs of K and of 
7 that are plausible, but arguably have a “solar 
system bias”. A systematic numerical study of the 
dynamical stability of multi-planet systems, over a 
wide range of planet masses and orbital periods is 
needed to improve the theoretically expected distri¬ 
butions of planetary system parameters and core¬ 
lations amongst them. Such a study requires sig¬ 
nificant computational effort, but is feasible with 


modern computers. This would enable an improved 
estimate of the planet mass function from observa¬ 
tional data of orbital periods alone, which are read¬ 
ily measured in almost all observational methods 
currently employed for the detection of extrasolar 
planets. 
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